function grav_init(params)


  % ****************************************************************
  %  output info
  % ****************************************************************

  % summary
  % disp(sprintf('Grid Size   : (%d, %d)', params.Nx, params.Ny));
  % disp(sprintf('Grid Length : %f', params.Lx));
  % disp(sprintf('Cell Length : %f', params.Dx));
  % disp(sprintf('Target Time : %f', params.Tfinal));
  % disp(sprintf('Time steps  : %d',  params.Tsteps));
  % disp(sprintf('Hydro Size  : %6.2fMB', 3*params.Nx.*params.Ny.*5.*params.Tsteps.*64/8/1024/1024));
  % disp(sprintf('Metric Size  : %6.2fMB', 3*params.Nx.*params.Ny.*5.*params.Tsteps.*64/8/1024/1024));



  % ****************************************************************
  %  prepare plotting
  % ****************************************************************

  if (params.plotMetric)
    % prepare figures for plotting
    close all;

    % load color maps
    load('cm', 'mycmap');
    global cmSize;
    cmSize = 512;

    % figures
    fig  = figure; 
    set(fig, 'Colormap', mycmap);
    set(fig, 'WindowStyle', 'docked');

    refresh();
    pause(0.001);
  end



  tic;

  Nx = params.Nx;
  Ny = params.Ny;
  Dx = params.Dx;
  Dy = params.Dy;

  % ****************************************************************
  %  compute derivative kernels
  % ****************************************************************
  % trefethen

  global d1x d1y d2x d2y;
  global gd1x gd1y gd2x gd2y;

  d1x = zeros(Nx, Nx);
  d1y = zeros(Ny, Ny);
  d2x = zeros(Nx, Nx);
  d2y = zeros(Ny, Ny);

  % x derivatives
  for x = 1:Nx
    for y=1:Nx
      s1 = 0;
      s2 = 0;
      for k = -Nx/2 : Nx/2
        s1 = s1 + 2. * pi * i / Nx / Nx * k * exp(2 * pi * i * k * (x - y) / Nx);
      end
      d1x(x,y) = real(s1)/Dx;
    end
  end

  % y derivatives
  for x = 1:Ny
    for y=1:Ny
      s1 = 0;
      s2 = 0;
      for k = -Ny/2 : Ny/2
        s1 = s1 + 2. * pi * i / Ny / Ny * k * exp(2 * pi * i * k * (x - y) / Ny);
      end
      d1y(x,y) = real(s1)/Dy;
    end
  end

  % compute chebyshev integration matrices
  global chebM1;
  global chebM0;

  N = params.Nz;
  D = cheb_init(N-1);

  ii  = 1:N-1;
  chebM1 = inv(D(ii,ii));

  ii  = 2:N;
  chebM0 = inv(D(ii,ii));

  disp(sprintf('Computed derivative kernels in %.2f sec', toc));
end
